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Abstract 

The transition to turbulence in Rayleigh-Benard convection with phase changes and the resulting 
convective patterns are studied in a three-dimensional Galerkin model. Our study is focused 
to the conditionally unstable regime of moist convection in which the stratification is stable for 
unsaturated air parcels and unstable for saturated parcels. We perform a comprehensive statistical 
analysis of the transition to convection that samples the dependence of attractors (or fixed points) in 
the phase space of the model on the dimensionless parameters. Conditionally unstable convection 
can be initiated either from a fully unsaturated linearly stable equilibrium or a fully saturated 
linearly unstable equilibrium. Highly localized moist convection can be found in steady state, 
in oscillating recharge-discharge regime or turbulent in dependence on the aspect ratio and the 
degree of stable stratification of the unsaturated air. Our phase space analysis predicts parameter 
ranges for which self-sustained convective regimes in the case of subcritical conditional instability 
can be observed. The observed regime transitions for moist convection bear some similarities to 
transitions to turbulence in simple shear flows. 
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I. INTRODUCTION 



Rayleigh-Benard (RB) convection, in which a fluid is heated from below and cooled from 
above, is one of the most comprehensively studied flow conflgurations, either in the regime 
of transition to convection and pattern formation 2| or in the fully turbulent case js, 4|. 
In nearly all these studies the working fluid is in a single phase. Moist convection combines 
Rayleigh-Benard convection and changes between two fluid phases, such as the gaseous and 
liquid ones. It is characterized by the internal release of latent heat for condensation and 
an increase of the efficient heat conductivity. This was demonstrated in a recent laboratory 
convection experiment with ethane [sj] and in numerical simulations joj. It is ubiquitous 
throughout the atmosphere of the Earth through the formation of clouds [2,8], but has also 
its importance in many industrial applications reaching from miniaturized heat exchangers 
9| to huge cooling towers in power plants jlo| . 

Condensation and evaporation are not only passive processes. The interaction of phase 
change and (turbulent) fluid flow results in a very complex dynamics intrinsic to the forma- 
tion of moist plumes and clouds. The thermodynamics of phase change is highly nonlinear, 
but can be considerably simplifled through the use of a piecewise linear approximation of the 
equation of state. This approach was suggested by Kuo 11 1 and Bretherton Q, Q and has 
been further developed recently jl^. Such simplifled thermodynamics enables a systematic 
way to investigate the dependence of convection on additional internal heating and to bridge 
the results directly to classical RB convection. This was done recently in direct numerical 
simulations (DNS) of turbulent convection with phase changes ISMlTj. However, even for 
a simplifled two-phase convection model the space of dimensionless system parameters is at 
least flve-dimensional and consequently prevents a detailed exploration by DNS. 

In this work, we present a systematic study of the dependence of moist Rayleigh-Benard 
convection on the dimensionless parameters. Our investigation is focused to the transition 
to moist convection. In order to perform such a study, we signiflcantly reduce the degrees 
of freedom in comparison to a fully resolved DNS. For this purpose, we construct a three- 
dimensional Galerkin model based on Fourier modes. This mode reduction limits our analysis 
to Rayleigh numbers of (9(10^) but makes it possible to monitor parameter dependencies 
systematically and to sample the structure of the underlying phase space by a statistical 
analysis in the spirit of recent investigations of the transition to turbulence in wall-bounded 
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shear flows 



18 



19|. 



Our focus is on moist convection in tlie so-called conditionally unstable (CU) regime. 
Dry air parcels are then stable with respect to vertical lifts while saturated air parcels are 
unstable and can rise in convective plumes. We address the following questions: How is the 
transition to convection in the CU regime characterized? How does it depend on the degree 
of saturation of the initial equilibrium state of the system? Which phase space structure is 
associated with the transition to moist convection? 

Models with reduced degrees of freedom are a powerful tool to study weakly nonlinear 
regimes in fluids. Several approaches to develop reduced models have been proposed, such 
as Galerkin models based on proper orthogonal decomposition 20|, reduced waveset approx- 



imations of turbulent flows 



23 



21 



22| or models in which only a least set of the largest-scale 



24| or a few dynamically relevant modes 25l-l27l| are captured. In dry Rayleigh-Benard 



convection such models provided a deeper understanding of the transition to convection 



23, l24] and the weakly nonlinear dynamics right above the onset of convective motion 
28|, l29|. Ogura and Phillips 30[ and later Shirer and Button 31[ applied systematic mode 
expansion to the moist convection case and set up minimal models. 

The transition to convective turbulence in dry Rayleigh-Benard convection proceeds via 
a sequence of bifurcations of increasing spatial and temporal complexity. Convection sets in 
with stationary roll patterns right above the critical threshold, bifurcates to other stationary 
states or to smooth time-dependent flow patterns for higher Rayleigh numbers and eventually 
becomes turbulent (see e.g. Ref. jl| for a comprehensive review). Dry convection is thus one 
of the standard applications of normal mode analysis to study linear stability and to seek 
for stationary nonlinear solutions at the transition threshold. Such normal mode analysis is 
only partly applicable for special conflgurations of moist convection since derivatives of the 
buoyancy fleld with respect to thermodynamic state variables have a discontinuity at the 
phase boundary, albeit the buoyancy fleld (which drives the fluid motion) itself is continuous. 

Bretherton analyzed the onset of moist convection by demonstrating the existence of 
linear two-dimensional growing modes and by studying their behavior in the weakly 
nonlinear regime of convection [13]. These studies showed that moist convection preferably 
develops in isolated saturated rising plumes separated from each other by broad unsaturated 
subsiding regions. This particular property of conditionally unstable convection was already 
predicted by Bjerknes [33] in a strikingly simple qualitative model. It is well-known from 
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field measurements such as in Refs. 



33 



34| and comprehensive large-eddy simulations 



of atmospheric moist convection containing additional physical processes such as radiative 



transfer, ice formation and precipitation 
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The purpose of our study is to extend Brethertons analysis of the conditionally unstable 
regime in various directions. It is based on the set of RB equations which we developed 
in 



jl^ and which are similar to [12|. First, we consider the three-dimensional case which 
has not been discussed in I2I, 1^. Second, the studies in 12, 13| are limited to a situation 



where the initial quiescent equilibrium is always exactly at the saturation line. This means 
that the layer is saturated but no condensed water is present in the initial profile, albeit 
subsequent convective motions can lead to further condensation. We will refer this regime 
to as the Kuo-Bretherton (KB) case. We extend conditionally unstable convection in both 
the subcritical case, corresponding to an initially unsaturated layer in which saturated air 
parcels may still be convective unstable, and the supercritical case, corresponding to an 
initially fully saturated layer in which condensed water is present at all levels This is 
achieved by increasing or decreasing the degree of stratification of the dry air compared to 
the parameters in the Kuo-Bretherton case. Third, we will perform a statistical analysis 
of the transition to turbulence that unravels the underlying phase space structure of our 
dynamical system in the sub- and supercritical cases. In particular, the subcritical case is 
known to be stable for small perturbations. A key issue here will be to determine for which 
ranges of the dimensionless parameters the self-sustained convective regimes can be obtained 
or not. We show that this transition - similar to the transition to turbulence in fundamental 
shear flows, such as a pipe or plane Couette flow - depends on the shape and the amplitude 
of the perturbation, a result of the competition between two co-existing attracting sets in 
phase space which are separated by a complex hyperplane. 

The outline of this paper is as follows. The model equations and the numerical implemen- 
tation of the Fourier- Galerkin approximation are explained in Section 2. Section 3 presents 
the two equilibrium states followed by section 4 that summarizes our analysis of the transi- 
tion behavior and the evolving convection states. Finally, we will give a short summary and 
outlook to future work. 
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II. GALERKIN MODEL 



A. Moist Rayleigh-Benard convection 



We now review briefly our model for moist 



Rayleigh-Benard convection. A detailed 



derivation and discussion can be found in Refs. |14j-|l6|. In the case of thermal convection, 
fluid motions are driven by a buoyancy force B that is added to the momentum balance 
equation to account for the variations of density in the fluid. The buoyancy field B is given 
by the following relation js^] 

B{S,qy,qi,p) = -g = , (1) 

with g being the gravity acceleration, p a mean density, p the pressure, S the entropy and 
qi the specific humidity of water vapor and liquid water. The following simplifications 
are made in the model. The first one is that the air parcels are in local thermodynamic 
equilibrium, which means that the two specific humidities can be combined to a specific 
humidity for the total water, a new state variable qx which remains constant for adiabatic 
non-precipitating processes: 

qT = qv + qi- (2) 

This reduces the number of prognostic variables to three. The second simplification is to 
apply the Boussinesq approximation to the system which in effect replaces the pressure 
dependence of 5 by a height dependence, i.e. B = B{S, qx, z). 

The dependence of the buoyancy field on the two remaining variables of state, S and qx, 
and the vertical coordinate z still contains the full thermodynamics of phase changes. One 
can approximate B{S, qx-, z) as a piecewise linear function of S and qr around the phase 
boundary between unsaturated gaseous and saturated liquid phases. This step preserves 
the discontinuity of partial derivatives of B{S, qr, z) with respect to qx and S at the phase 
boundary (and thus the release of latent heat). The last two steps limit the applicability 
of the model to shallow convection where the thickness of the atmospheric layer remains 
moderate. 

Since i? is a linear function of the two variables of state S and gr, we can introduce 
two new prognostic fields, a dry buoyancy field D and a moist buoyancy field M as linear 
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combinations of S and qt- An air parcel is saturated and liquid water is present when [l^ 

M -D + N^z>0. (3) 

Ns is the Brunt- Vaisala frequency which is defined as = g(Td — Ts)/Tref- Here, and 
r<j are the dry and moist adiabatic lapse rates and T^e/ is a reference temperature, e.g. the 
temperature at the bottom plane. All these assumptions lead to a formula of the buoyancy 
field B{M, D, z) which is given by 

5(x, t) = max (M(x, t), D{x, t) - N^z) . (4) 

The saturation condition (jl]) is a nonlinear relation which can be easily carried out in nu- 
merical simulations. On the basis of (jl]), we can define the amount of condensate by 

gz(x, t) = M(x, t) - D(x, t) + Nlz > . (5) 

Consequently, all points > belong to moist plumes with droplets or clouds. 
The dry and moist buoyancy fields can be decomposed into 

L'(x,t) = 'D{z) +D'{^,t) 

= Do + ^^^^^z + D'i^,t) (6) 

M(x,t) = M(2) + M'(x,t) 

= M, + ^^^^z + M'i.,t). (7) 

£1 

Here, Dq and Mq are prescribed values at 2; = respectively, Dh and ai z = H. The 
variations about the linear equilibrium profiles of both fields, D' and M', have to vanish at 
z = and H. Equation (jl]) can now be transformed into 

B = M{z) + max (M', D' + D{z) - M{z) - N^z) . (8) 

To obtain a dimensionless version of the model equations, one has to define characteristic 
scales. These are the height of the layer H, the velocity Uf = a/ H{Mo — Mh), the time 
T = H/Uf, the characteristic kinematic pressure Up and the buoyancy difference Mq — Mh- 
The dimensionless equations contain three nondimensional parameters, the Prandtl number 
Pr = u/k and the dry and moist Rayleigh numbers Rao and RaM which are given by 

UK UK 



Here v is the kinematic viscosity and n is the diffusivity of the buoyancy fields. In addition, 
two more parameters arise from the saturation condition ([8]): the Surface Saturation Deficit 
{SSD) and the Condensation in Saturated Ascent (CSA) 

SSD= CSA=JSIL^. (10) 

The parameter SSD determines the degree of saturation of the air parcels at the bottom 
plane z = 0. For the rest of the work, we set Dq = Mq and thus SSD = 0. The second 
new parameter, CSA, is proportional to the amount of water formed in a saturated parcel 
adiabatically rising from the bottom to the top. The buoyancy B in dimensionless form is 
given by 

B = max (^M', D' + SSD + - - CSA^ . (11) 

We can decompose the buoyancy field in the last equation in a mean contribution B{z) which 
can be added to the kinematic pressure p and a fiuctuating field B', i.e. d^p + B = d-^p + B' 
. The dimensionless equations are then 



dtu + (u • V)u = -Vp + ^T&V^u + B'e, (12) 

V ■ u = (13) 

d,D' + (u . V)D' = ^^V^D' + f^n. (14) 

d,M' + (u ■ V)M' = -^=2=V'M' + uz . (15) 

Here u is the velocity field. Together with ( [TT]) this forms a closed system of equations 
that describes moist Rayleigh-Benard convection. Beside SSD = 0, we set Pr = 0.7 and 
CSA = 4/3 in this work. 

An alternative choice for the two non-dimensional parameters which are related to the 
phase change is given by 

SSD 



CWo 



CWh = I + CW0 + 



CSA 

Ran/RaM — 1 



CSA 



A positive value of either CWq or CWh corresponds to cloud water (CW) present at the 
bottom or top. A negative value stands for a water deficit. Here CWq = 0, while CWh will 
vary with a change of Rao- 
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B. Galerkin approximation 



The model equations are solved in a rectangular box of height 1 and width AR, which 
is the aspect ratio of the domain. In the horizontal directions x and y we use periodic 
boundary conditions. In the vertical direction z we apply free-slip boundary conditions at 
z = 0,l 

u, = D' = M' = 0, ^ = 1^ = 0. (16) 

Galerkin approximations are frequently used to investigate dynamics that mainly includes 
a few modes 
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27 



31 



381]. We will use this method here for comprehensive statistical and 
parametric studies that would not be possible with direct numerical simulations in which 
all modes to the smallest flow scale are included. Our Galerkin approximation is based on a 
truncated Fourier series. For a simple implementation of the boundary conditions we double 
the box from z = —1 to z = 1. For a wavevector n = (n^;, n„, ra^) G Z'^ we introduce the 



vrm^ . (17) 



following notation 

" ' AR ' AR 

Note that the aspect ratio is included in this definition and does not emerge directly in the 
following equations. Every field is expanded in a Fourier series. For example the velocity 
field is given by 

u (x, t) = ^ u (n, t) e"'' . (18) 
The following set of ordinary differential equations results from f[T^ - 0151] 
dtuin.t) = - J2 (u(p,t)-q)u(q,t) + 5'(n,t)e, 

p+q=n 



Pr 

-np(n,t) + W^nVn,t) (19) 
u(n,t) ■ n = (20) 

yjFrRaM Ro-m 

- J2 (u(p,t)-q)/^'(q,t) (21) 

p+q=n 

dtM'{n,t) = M'{n,t)+u,in,t) 
^FrRau 

- J2 (u(p,t)-q)M'(q,t). (22) 

p+q=n 
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Since the number of modes is small compared to a direct numerical simulation, the convolu- 
tion sums in Eqns. f|T9|) - fl22|) can be calculated directly. Taking the divergence of ([T9l) and 
using (I2CT]) we get the following equation for the pressure, 



np(n,t) = ((u(p, t) • q)(u(q, t) • n) 



p+q=n 

-ni7m^B'{n,t)) . (23) 



Since all quantities have to be real in physical space it follows that 



u(-n,t) = u*(n,t) (24) 
M'(-n,t) = M'*{n,t) (25) 
D'{-n,t) = D'*{n,t), (26) 



where the asterisk stands for complex conjugate. The free-slip boundary conditions (fT6|) 
are implemented by an additional symmetry with respect to the plane z = for u^, Uy and 
antisymmetry for m^, M' and D' . This results to 

Ux{x,y,-z,t) = Ux{x,y,z,t) (27) 
Uz{x,y,-z,t) = -u^{x,y,z,t) . (28) 

In Fourier space this requires 

u^{n^,ny,-n^,t) = ul{n^,ny,n^,t) (29) 
u^{n^,ny,-nz,t) = -ul{n^,ny,n^,t) . (30) 

Bretherton [12I] noticed already that a further simplification can be made by using that 
M' and D' obey the same linear equations (except for a different driving amplitude) and 
boundary conditions. We can consequently use 

D' = ^M'. (31) 

Thus only one buoyancy field has to be simulated. The Fourier coefficients B'{n,t) are 
calculated in the following way. First B' is evaluated in physical space on a grid. Then the 
coefficients 

5'(n,t) = :i /" 5'(x,t)e""rfx (32) 
^ Jv 
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FIG. 1: Examples for runs in the subcritical (solid lines) and the supercritical cases (dotted lines). 
In each regime we perturb the equilibrium either with a finite initial perturbation or an infinitesimal 
perturbation. The transition to convection or the return to the equilibrium are displayed by the 
corresponding kinetic energy versus time. 



are computed using 4th-order rule for the evaluation of integrals [39|. The time stepping 
is done with a third order Runge-Kutta scheme by von Heun [4^. The choice of the wave 
vectors is a compromise between computational costs and incorporating the large-scale flow 
aspects. Note also that the saturation condition (ITT]) will generate a fair amount of scale 
interactions. Furthermore remember that B' is not differentiable at the phase boundary, but 
approximated by continuous differentiable functions. Therefore a certain number of modes 
is necessary to get an appropriate approximation. All Galerkin model runs at AR = 3,4 
are done with the wave vectors {\nx\ < 5, \ny\ < 5, \nz\ < 5) if not stated otherwise. Mode 
(0,0,0) is always excluded. For larger aspect ratios more wavevectors are taken. These are in 
total 975 independent degrees of freedom since one velocity component can be determined 
by the incompressibility condition. 

III. CONDITIONALLY UNSTABLE EQUILIBRIA 

In the absence of fluid motion, the profiles of D and M will be linear across the layer 
as given in Eqns. (jS]) and ([7]). The choice of the four amplitudes, Dq, Dh, Mq, and Mh, 
causes different classes of equilibrium states. Conditionally unstable layers are defined as 
layers in which clear air is stably stratified, i.e. > Dq and thus Ra£) < 0. Moist air 
is unstably stratified, i.e. < Mq and thus Rgm > 0. To classify the CU equilibria of 
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a moist convection layer one defines the Convectively Available Potential Energy (CAPE) 



ill by 

CAPE= [ [B(D{0),M{0),z) - B(D{z),M{z),z)]dz. (33) 
Jo 

The first term describes the potential energy of the air parcels that start at the bottom. 
The second term stands for the potential energy of the background equilibrium profile. A 
necessary condition for the onset of convection is that CAPE > 0. This is equivalent to 
{SSD = 0) 

CAPE > ^ CSA > (34) 

In the following, we will discuss two cases of conditionally unstable convection and investigate 
their stability properties. 

The quiescent equilibrium configuration of the convection layer is fully subsaturated if 
the linear profiles satisfy 

M{z) < D{z) - N'^z , (35) 
for all z E [0, 1]. This results together with ([TTj) to 

CSA>-^>CSA-1. (36) 

We will refer convection which arises from this equilibrium case to as the subcritical CU 
regime. In this regime, the linear profile is unsaturated, and small perturbations cannot 
change the saturation of air parcels. In effect, small perturbations experience the stable 
stratification of the dry buoyancy field D, and the diffusive equilibrium solution is linearly 
stable. Finite perturbations however can trigger moist convection out of the equilibrium 
(see next section). This is demonstrated in Fig. [T] which compares the evolution of the 
kinetic energy for two perturbations of different initial amplitudes. Note that if the stable 
stratification of the dry buoyancy field is too strong, CAPE cannot be positive. This is 
the absolute stability threshold, a necessary (but not sufficient) condition to initiate moist 
convection. It defines one interval boundary of the subcritical regime. 

The other interval boundary of the subcritical regime is set by the so-called Kuo- 



Bretherton (KB) case Uj, ll2|. This case is established if the degree of stable stratification 



of D is reduced until M{z) = D{z) — N^z corresponding to 

- Rao/Rau = CSA - 1 . (37) 
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The whole domain is held in an equilibrium at the saturation threshold. This case has been 
extensively studied by Bretherton [12I, 13 1 who demonstrated the existence of exponentially 
growing solutions in the form of isolated saturated plumes and investigated their behavior 
in the weakly nonlinear regime. 

Once the stratification is further reduced, small-amplitude perturbations to the equilib- 
rium can trigger moist convection as demonstrated by the dotted curves in Fig. [1] The 
quiescent equilibrium is fully saturated for 

M{z) > D{z) - Nlz . (38) 

We will refer convection triggered from this equilibrium state to as the supercritical CU 
regime because the equilibrium is linearly unstable, which corresponds to non-dimensional 
parameters satisfying the condition 

CSA-l>-^>0. (39) 

Note that when the value of CSA is less than 1, there is no supercritical CU regime due to 
the fact that saturation of the entire layer would require a positive value of the dry Rayleigh 
number Rao for these values of the moist Rayleigh number Rqm and CSA. The case of 
both buoyancy fields unstably stratified, corresponding to Ra^ > and RaM > 0, is also 



possible and has been discussed in detail in 
both CU regimes. 



14j-|l6|. In the following section we will discuss 



IV. TRANSITION TO MOIST CONVECTION 



A. Subcritical regime 



1. Statistical analysis of the transition to convection 

In the subcritical regime, the linear equilibrium solution is stable. After an initial per- 
turbation and burst of convection, the layer can either return to the diffusive equilibrium 
state (denoted as A in the following) or ends up in a convective state with fluid motion 
(denoted as B) for the same set of dimensionless parameters RaM, Ran and AR when a 
finite perturbation of different shape or amplitude is applied. Such a behavior is similar to 
the transition to turbulence in wall-bounded shear flows p^]. To investigate this behavior 
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FIG. 2: Statistical analysis of the transition to moist convection. We show the transition proba- 
bilities to obtain a convective solution B. Left: Dependence on RaB for fixed magnitude of initial 
perturbation, E]^in fmd E^' ■ The vertical line marks the dry Rayleigh number corresponding with 
the KB case. Right: Amplitude dependence of the transition probability. Solid line is for transi- 
tion from convective state B to diffusive state A, dotted line from state A to state B. The vertical 
line marks the kinetic energy of state B for i^au = -1.60 x 10"^. AU data are for AR = 4 and 
Ro-M = 3.73 X 10^ which is for the subcritical CU regime. 



systematically, we fixed the aspect ratio to 4 and performed a statistical investigation in the 



spirit of those in shear flows 
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431]. Therefore, we took 192 randomly selected perturba- 



tions among the wave vectors with wave numbers In^l = 1, 2, \ny\ = 1, 2, and jn^l = 1,2 for 
the velocity field and the moist buoyancy field. The initial kinetic energy E^in = (u^)v/2 
and the variance of the moist buoyancy field Em' = {M''^)v/2 were the same for all 192 
runs. The result is shown in Fig. |2](left). For Rao chosen in the supercritical regime, state 
B is always obtained. The probability to end in steady moist convection is unity. For RaD 
taking values in the subcritical regime, the probability of a relaxation to state B decreases to 
zero with decreasing Ra£,. Note that the value of Rac with zero probability to obtain B is 
still significantly larger then the dry Rayleigh number which is associated with the threshold 
of CAPE = 0. This would result to Rao = —4.98 x 10^ for the present example. 

We also determined the smallest value of Ra^ to observe a (stationary) convection state 
B. State B for Rao close to zero transition probability is taken as an initial condition for 
a new simulation with a slightly reduced Ra£,. This case is advanced in time as long as it 
needs to relax anew. The procedure can be continued iteratively and the relaxation to the 
new (stationary) state increases in time for decreasing Rao. In this way we found state B 
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FIG. 3: Transition to moist convection as a function of the perturbation amplitude and the dry 
Rayleigh number. Left: AR = 4. Here, 120 different values for Rao and 192 different amplitudes 
are taken. Right: AR = 3. Here, 40 different values of Rao and 94 different perturbation 
amplitudes are taken. Gray stands for reaching state B out of the equilibrium, white for a return 
to state A. The vertical line marks the KB case at the border between sub- and supercritical 
regimes. Each simulation was run for 300 dimensionless time units which is more than a diffusive 
time td = H'^ /v. 



appearing at Rayleigh numbers as small as Ran = —1.99 x 10^. 



2. Dependence on initial perturbation amplitude 

Results on the amplitude dependence of the transition at a fixed set of system parameters 
are shown in Fig. |2] (right). We started therefore either from state A or B and determined 
the probability to reach state B or A, respectively, by varying the amplitude of the initial 
perturbation. The moist buoyancy field M' is not perturbed. One can observe that the 
convective state B is very stable with respect to small perturbations (solid line) and becomes 
unstable when the kinetic energy of the perturbation exceeds a threshold which is marked 
as a solid vertical line in the figure. In contrast, the diffusive state A (dotted line) behaves 
differently and shows a finite probability to switch to the convective state B over nearly the 
whole range of perturbations. 

The observation suggests a more systematic monitoring of the plane which is spanned 
by the perturbation amplitude and RaD, as done for example in a Galerkin model for a 



plane shear flow 



271]. Now the shape of the perturbation is fixed. The amplitude of the 
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perturbation and the dry Rayleigh number are varied to cover both CU regimes. The detailed 
scan of the parameter plane is shown in Fig. |3] for aspect ratios of 3 (right) and 4 (left). 
Perturbations that cause a return to A are in white and those that initiate a convection 
state are in gray. 

The figure illustrates nicely that the boundary between the two attracting regions is 
complex if one samples the phase space, even for the moderate Rayleigh numbers discussed 
here. This is for example obvious from the isolated island that can be seen in the left panel of 
Fig. [3] for AR = 4. The observation is similar to what has been observed in low- dimensional 
shear flow models. An additional property of shear flows is the finite lifetime of transient 



states 
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43| which we did not observe for the present system. Once moist convection is 



initiated it remains sustained, except at the boundary of the basin of attraction of state B 
(see Fig. |6]and subsequent text). 



3. Structure of convection states in the subcritical regime 

When convection is initiated, state B takes the shape of a single localized steady moist 
plume (or cloud) at small aspect ratio. It is characterized by a strong saturated upward 
motion inside the cloud balanced by a weak downward motion outside - a well known feature 
of moist convection in this regime as already mentioned in the introduction and e.g. in 
The position of the single cloud in the box differs from case to case, but the shape and the 
flow structure are always the same. It is observed for moderate magnitudes of Ra£, and 
Rqm at different aspect ratios as seen in Fig. HI Similar solutions have been reported by 
Bretherton 
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13| in his two-dimensional simulations of the nonlinear evolution. 
For AR = 2, the convective state B was not detected at all. The size of the box is then 
too small to provide enough space for the dry subsidence outside the cloud and thus to form 
a stable cloud pattern. For the remaining cases, i.e. AR = 3, 4, 5, 6, a single cloud appears 
in the box. Neighboring clouds are stabilizing each other (recall the periodic boundary 
conditions in x and y). States B for several AR do not differ qualitatively in their overall 
structure as supported by the mid and bottom panels of Fig. |H With increasing AR, the 
clouds become more and more smooth and grow in diameter relatively to the aspect ratio. 
Quantitatively, the velocity magnitude and the cloud fraction (all points with qi > 0) of the 
layer increase with the aspect ratio, as displayed in Fig. |5l 
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FIG. 4: Examples for stationary single moist plume (or cloud) at Rao = —1.50 x 10^ and Rum = 
3.73 X lO'*. The upper plot shows the kinetic energy vs. time that becomes eventually stationary. 
The solid line corresponds to AR=6 and the mid plot, the dashed one to AR=4: an the lower plot. 
The mid and bottom panels show centered cross section contour plots in the x — z plane of the 
final (stationary) states. We show the velocity field vectors in the plane and the contours of liquid 
water content qi as defined in dSJ. The bright solid line in the mid and bottom figures marks the 
boundary with qi=0. All points inside this boundary with qi > belong to the cloud. Points with 
qi < are outside the cloud and stand for a water deficit. 
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FIG. 5: Maximum vertical velocity inside the moist plume and fraction of saturated air with qi >0 
as a function of the aspect ratio. Here, Rao = —1.50 x 10^ and RaM = 3.73 x 10^. 

For even larger aspect ratios, such as AR = 7 and 8, we still observe a single cloud in the 
domain, but now temporal fluctuations of all fields can evolve. The stabilizing mechanism 
from neighboring clouds is too weak. It is mentioned here that the number of wave vectors 
and the computational grid have been always adjusted with growing AR in the Galerkin 
model. Multiple clouds with periodic distance apart have also been observed in 

4- Equilibrium and moist plume as coexisting attractors 

One possible explanation for the observed transition behavior in the subcritical regime is 
that two coexisting attractors (or fixed points) which are associated with equilibrium state 
A and (stationary) convective solution B. Both attractors and their associated basins of 
attraction in phase space will vary with changes in the parameters. While the basin of B 
will grow with increasing Ran, the basin of state A will shrink. The coexistence of two 
attractors is also supported by the observation that initiated convection states do not decay 
after a finite time. 

To confirm that there is a small attracting set in the phase space close to both, states A 
and B, we proceed as follows for one particular parameter set. Both fixed point solutions, 
denoted as and X^, are slightly perturbed. The time advancement is started with 
initial conditions taken from a sphere in phase space around both attracting states, i.e., 
X-^ ^(t = 0) = X^ Q + ^X-' with j = 1, 1016. They are randomly selected on this sphere 
with I^X-'I = 5 X 10~^. All cases relaxed either to states A or B. The dry Rayleigh number 
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was fixed at Rao = —1.85 x 10"^ for this investigation, and the aspect ratio is AR = 4. 
Our study is not a rigorous proof, but provides evidence for our present picture of phase 
space. Repeating this procedure in the supercritical regime at Ra^ = —10^, we find that the 
system escapes from state A for all 1016 initial conditions, while state B remains attractive 
as before. In other words, the saturated equilibrium state (state A) is now linearly unstable. 



5. Tracking the edge state between both attractors 

In the following, we explore the boundary between both basins of attraction in more 
detail. Solutions which can be found on this boundary are denoted as edge states [4J, |45 |. 



We used an edge tracking method which is similar to Ref . [46| . The boundary that separates 
both basins of attraction must be somewhere between states A and B, and can be obtained 
ClS db linear combination of the two fixed points: 

X{t = 0) = XX*^ + (1 - X)X*B , (40) 

with A G [0,1]. Because state A is the origin it reduces simply to AX^. We divided 
the interval [0,1] into 65 equidistant points and detected the bin with the value of A that 
separates state A from state B. This procedure is refined successively and the refined values 
of A are used to start a new run at t = 0. Results for Ra^ = —1.60 x 10^ are shown in Fig. 
[61 The detected edge state is a periodic orbit. It is a cloud fixed at one position that first 
grows to a cloud with shape and flow structure very similar to state B (see Figs. H] and [8] 
(top)). Afterwards it collapses to a small cloud pen, but never vanishes. 

A picture of the phase space structure in the subcritical case is illustrated in Fig. [7l Shown 
are the original equilibrium state A, the steady moist convection state B and the edge state. 
The figure indicates also the complex interface between both basins of attraction that we 
detected in several ways in the present study. 



B. Supercritical regime 



1. Structure of convection states in the supercritical regime 

In the supercritical regime, the diffusive equilibrium is saturated and linearly unstable. A 
convective state can always be initiated. This was already shown in Fig. [3] for aspect ratios 
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FIG. 6: Edge state tracking for AR = 4 and RaD = —1.60 x 10^. We plot the kinetic energy versus 
time. Shown are always the two runs at the given refinement level that enclose the edge state. 
To distinguish the refinement levels the corresponding kinetic energy curves are plotted dotted 
and dashed in alternating sequence. The last iteration that was possible with double precession 
corresponds to the solid curves. 




3 (right) and 4 (left). Typical convection states B in the sub- and supercritical regime are 
compared in Fig. |8l One can see that the shape of the moist plume changes qualitatively. 
While in the subcritical regime an isolated cloud exists, in the supercritical regime a closed 
cloud layer at the top is found. For decreased stratification the moist updraft becomes 
stronger while the overall shape remains the same. In Fig. |9] (top), the L2 norms of the 
velocity components as a function of Rqd are shown. The two horizontal velocity magnitudes 
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FIG. 8: Velocity and liquid water content in the x — z plane for AR=4. The bright line marks again 
qi = 0. Top figure: Rajj = — 1.96x 10^ in the subcritical regime. Bottom figure: Rao = — 7.85x 10^ 
in the supercritical regime. For both cases Rum = 3.73 x 10^. 

are the same as one would expect due to the symmetry. All three contributions to the kinetic 
energy grow steadily until state B becomes non- stationary in the supercritical regime. This 
is observed at about Rao = —7820 (see Fig. [TOj) . Interestingly, the maximal vertical velocity 
inside the cloud varies much less for the range of Ra^,. The bottom panel of the same figure 
displays the L2 norms of moist and dry buoyancy. The moist up draft increases in intensity 
as we reduce the stable stratification which is indicated by ||M'||2. The shrinking of the 
buoyancy fluctuation at the KB case is due to the fact, that the background equilibrium 
switches from fully unsaturated to fully saturated case. Note also that the moist Rayleigh 
number Rqm remains constant in this study. 

2. Recharge- discharge convection 

For dry Rayleigh numbers in the supercritical regime right above the KB case, we observe 
in some instances moist convection in a highly intermittent recharge-discharge regime which 
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FIG. 9: Dependence of several fields on RaD for convection in steady state B in the cell with 
AR = 4. The dotted vertical line marks the KB case right between subcritical and supercritical 
regimes. We take the L2-norm of the fields which is given by || • ||2= l/V Jy | • pc?x. For the whole 
range of Rao convection is still in a stationary regime. 
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FIG. 10: Transition from stationary to time-dependent moist convection. The dry Rayleigh number 
Ran is therefore increased from -7840 up to -7540 in increments of 60. Data are for AR = 4, 

SSD = 0, CSA = 4/3, and Rum = 3.73 x 10^. 
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FIG. 11: Example for the recharge-discharge convection for AR = 4. The initial equilibrium state 
is fully saturated and very close to the KB case. When time axis is rescaled by the diffusion time 
of the flow, td = H^/j^, the two recharging events will appear at t/td = 1 and 2. 

is displayed in Fig. [11] via the graphs of the L2 norms of the velocity and the buoyancies 
versus time. After an initial convective burst of activity, the system slowly relaxes toward 
the quiescent equilibrium. However, after a while, convection is reinitiated and decays 
again. Interestingly, a similar behavior has been found in a channel flow with a strong 
transverse magnetic field that enforces a switching between two-dimensional and three- 
dimensional flow states. The reason of why the conditionally unstable moist convection is 
not fully relaxing to the quiescent state can be explained as follows. In the supercritical 
regime, the equilibrium state is unstably stratified. Small perturbations can trigger new 
moist convection. We will get back to this discussion at the end of this section. 

The number of such cycles can vary from one to several. An observed trend is that 
the closer the parameter set is to the KB parameter set the more cycles are possible. In 
the Galerkin model, it is found that the recharge-discharge regime also depends on the two 
Rayleigh numbers. For higher Rayleigh numbers, more cycles can be observed, which implies 
that viscosity plays an important role for this phenomenon. In particular, the duration of the 
cycle is approximately equal to the diffusive time-scale td = W^/v. This indicates that the 
quiescent relaxation is associated with the diffusion of dry and stably stratified air across 
the entire layer. However, as the origin is a linearly unstable equilibrium, the quiescent 
relaxation cannot be sustained indefinitely. It is terminated by the abrupt onset of a cloudy 
moist plume. This behavior has been reproduced in direct numerical simulations (DNS) 
which will be reported elsewhere jirj l. 
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FIG. 12: Snapshots of the hquid water content qi and the projection of the velocity field into the 
plane for the same data as in Fig. [TTJ They are taken in the plane with angle 7r/4 to the x — z 
and y — z plane. Bright line marks qi = 0. From top to bottom, first panel is at t = 183. Velocity 
vectors are enhanced by a factor of 250. Second panel is at t = 198 and velocity enhanced again 
by 250. Third panel is taken at t = 210. The velocity vectors are now stretched by a factor of 6. 
Fourth panel taken at t = 224, velocity is enhanced by a factor of 4. Note that the colorbars of 
the panels differ by orders of magnitude. 
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We illustrate this dynamics in detail for an example in Figs. [TT] and [121 Different phases of 
the evolution are displayed. The cycle starts at t = 27 and is characterized by an exponential 
decay of all quantities (see Fig. [TT]) . At that time, the entire box is filled with unsaturated 
air and is stably stratified, except in the vicinity of the upper plate which is weakly saturated 
due to the boundary conditions. Since the system converges to the quiescent equilibrium 
state, moist air starts filling slowly the box from the upper plate. This process begins at 
t = 171 where the buoyancy B' is not decreasing any further because of the occurrence 
of moist buoyant air. The first panel of Fig. [12] taken at t = 183 illustrates this stage 
slightly before the onset of convection. Fluid motions are very weak, and mostly horizontal. 
Saturated air begins to fill the layer slowly from above. At time t = 191, the buoyancy is 
increasing rapidly and the moist air in the box is rearranged to column-like clouds as shown 
in the second panel for t = 198. The solid bright line is again the cloud boundary. The 
velocity remains small. But at t = 203 the buoyancy is strong enough to form moist upward 
motion in the cloud columns (see the third panel at t = 210). Due to latent heat release 
this process is self-amplifying. The velocity increases very rapidly and the cloud columns 
form intense single updrafts up to t = 224 (fourth panel of Fig. [12]) . If the updrafts are 
too strong, the convection can over-stabilize the layer which then must be destabilized by 
diffusion again. This means that a new recharge-discharge cycle is initiated. Alternatively, 
weaker updrafts may lead to a sustained convective regime in which the stabilization by 
convection balances diffusion on the same time-scale. In general, the larger the domain the 
easier to escape this cyclic recharge-discharge convection. 

In order to detect which Galerkin modes are mainly included in the recharge-discharge 
dynamics, we reduced the wavenumber space stepwise. It was possible to reproduce the 
cycles with wavenumbers (|?T,a;| < 1, < 1, I^t-^I < 2) (see Fig. [T3]) . We can thus conclude 
that this regime of convection is dominated by the large-scale degrees of freedom in the flow 
model. 

The conditions for the onset of a recharge-discharge cycle can be understood as follows 
(see also Fig. [12]) . First, nearly the whole box must be filled with dry and stable air. Second, 
the velocity field is then almost horizontal, because an upward motion could create moist 
and unstable air and thus more upward motion. If we take such an initial condition for less 
stably stratified dry air (i.e. a lower Rao) we can produce at least one cycle. For parameter 
sets close to the KB case it is easier to satisfy these conditions. A small perturbation in 
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FIG. 13: Effect of mode reduction on the recharge-discharge convection regime. Solid line is for 
the usual resolution in our Galerkin approximation, dashed line for the reduced resolution. 



FIG. 14: Sketch of the phase space structure which is associated with the recharge-discharge 
behavior. The gray arrow marks the case of a perturbation only in the Ux and Uy direction. The 
dashed line is a example for a two cycle recharge-discharge run. 

the dry buoyancy field is enough to get an almost dry domain. That is why we can observe 
this phenomenon mostly close to the KB case. The closer we are to the KB case the more 
cycles are observable. In Fig. [H] we give a picture for the phase space structure for the 
recharge-discharge convection. It is clear from the equations of motion that perturbations 
in Ux and Uy only will always relax to state A. This path corresponds with a stable manifold 
and is indicated by the gray arrow in Fig. [HI 
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V. SUMMARY AND DISCUSSION 



The transition to turbulence in moist Rayleigh-Benard convection starting from a con- 
ditionally unstable (CU) equilibrium has been studied in a Fourier- Galerkin model. We 
classified the initial CU equilibrium of the layer in two different regimes, the subcritical 
unsaturated and the supercritical saturated. The transition behavior and the evolving dy- 
namics are significantly altered by the degree of saturation which is controlled at the top of 
the layer by Rao (when the parameters SSD, CSA, RaM and Pr are held fixed) and the 



aspect ratio AR of the layer. Our study extends earlier results of Bretherton [1^, Il3| by 
demonstrating the existence of nonlinear three-dimensional convective regimes in a broad 
portion of the parameter space. 

In contrast to classical dry RB case, convection is localized and does not fill the entire 
layer with rolls and thermal plumes. In the subcritical case, we showed that a narrow moist 
plume is surrounded by dry unsaturated air in which the fiuid moves downward by diffusion 
(see Fig. Hj). This characteristic structure is well-known from observations and numerical 



simulations 132, 
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ll|-|l3l|. In the supercritical 



I - also from those of the KB regime 
regime, the top layer is however always saturated, a closed cloud layer is then fed by a 
narrow moist plume from below (see Fig. [8]). The up-down symmetry as known from the 
dry RB case is consequently broken. Moist CU convection is additionally constrained by 
the aspect ratio. Time- dependent and eventually turbulent convection requires sufficiently 
extended cells. 

In the subcritical case, when the atmosphere remains subsaturated at the top {M{z) < 
D{z) — NgZ), two stable attractors or fixed points coexist. The first one corresponds to a 
diffusive atmosphere with no motion, while the second one exhibits an overturning circu- 
lation in the form of a steady rising moist plume (or cloud) balanced by subsidence in the 
unsaturated environment. The existence of these multiple attractors has been confirmed 
through several statistical analyses. With decreasing Rao, it becomes less likely that a con- 
vective equilibrium can be obtained. For a sufficiently negative dry Rayleigh number, that is 
however still well above the threshold for which the layer is absolutely stable [CAPE < 0), 
the probability to switch to moist convection becomes zero. In this range of Ra^, we traced 
a periodic orbit at the edge of the basin of attraction of the convective state (see Fig. [7j). 



Some of the transition properties in this regime are similar to those of wall-bounded shear 
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flows. We mention the sensitive dependence on the amphtude of the finite perturbation 
and the complex shaped boundary between laminar and turbulent states. Furthermore the 
convective state is localized in space, very similar to turbulent spots in a shear flow 48 |. 
But there are also differences. In the present model, we did not observe transient convection 
states which suggests the existence of an attractor rather than a chaotic saddle (see Fig. [7j). 
An exception is the edge state. 

If the layer is held saturated at the top {M{z) > D{z) — N^z), the initial equilibrium is 
linearly unstable with respect to small amplitude perturbations. The system switches into 
a moist convection mode which in dependence on the chosen aspect ratio and the degree of 
saturation at the top can be a recharge-discharge mode or a stationary moist plume mode. 
A further increase in the dry Rayleigh number Ra^i, which corresponds to a decreasing 
stratification in the unsaturated environment, causes the moist plume structure to become 
time-dependent and eventually turbulent. 

A few more words are in order now. The present analysis is focused on moderate 
Rayleigh number and small aspect ratios of 3 to 8. These choices, made necessary by 
the Galerkin truncation, will have their influence on the initialization of cloud patterns. 
Moderate Rayleigh numbers and thus larger viscosity and diffusivity ease the onset of sus- 
tained convective motion outside the cloud since the stable dry air can descent faster there 
by diffusion. With increasing Rayleigh number this process becomes increasingly inefficient. 
In other words, the upward transport in ascending plumes must be balanced by slower sub- 
sidence occupying larger and larger portions of the domain around the clouds. Thus, bigger 
aspect ratios are likely to be necessary to obtain a self-sustained convective regime at high 
Rayleigh number. This is confirmed by recent high-resolution DNS of this model at higher 



Rayleigh numbers 



ITj . Note also, that with a view to moist convection in the atmosphere. 



radiative cooling or non-equilibrium effects such as precipitation step in and can become im- 
portant physical processes to establish moist convection in the conditionally unstable regime. 
They provide the additional dissipation mechanisms that amplify the downward fluid trans- 
port outside the moist plume when diffusive transport becomes increasingly inefficient with 
growing Rayleigh (or Reynolds) numbers and will ease the constraint of larger aspect ratios. 

The results presented here can be expanded in several directions. First, our study only 
considered the case where the lower boundary is exactly at the saturation boundary, cor- 
responding to SSD = 0, which can be interpreted as convection over an ocean surface. 
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Further exploration of state space to include cases with non-zero SSD could shed some 
light on the different behavior between convection over land and over the oceans. Second, 
changing boundary conditions from free-slip to constant flux conditions can have an^infiuence 
on the dynamics at moderate Rayleigh number as known from dry convection [49]. Third, 
the addition of a radiative cooling, which occurs in the atmosphere due to the emission of 
infrared radiation by greenhouse gases, will significantly enhance the number of dynamical 
regimes as has been recently discussed by large eddy simulations with a parametrization of 
small-scale turbulence in Ref. [50 1. 
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